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Abstract 

Considering the coherent nonlinear dynamics in double square well potential we find the example 
of coexistence of Josephson oscillations with a self-trapping regime. This macroscopic bistability 
is explained by proving analytically the simultaneous existence of symmetric, antisymmetric and 
asymmetric stationary solutions of the associated Gross-Pitaevskii equation. The effect is illus- 
trated and confirmed by numerical simulations. This property allows to make suggestions on pos- 
sible experiments using Bose-Einstein condensates in engineered optical lattices or weakly coupled 
optical waveguide arrays. 

PACS numbers: 05.45.-a, 42.65.Wi, 03.75.Lm 
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I. INTRODUCTION 



The problem of nonlinear dynamics in double well potential has been first addressed by 
Jensen^ who considered light power spatial oscillations in two coupled nonlinear waveguides, 
which resemble Josephson oscillations^ in the spatial domain. The latter macroscopic quan- 
tum tunneling effect, originally discovered in superconducting junctions, is caused by the 
global phase coherence between electrons in the different layers. More recently the similar 
realization of a bosonic Josephson junction has been reported for a Bose-Einstein conden- 
sate embedded in a macroscopic double harmonic well potential^. The difference with the 
ordinary Josephson junction behavior is that the oscillations of atomic population imbalance 
are suppressed for high imbalance values and a self-trapping regime emerges^. 

The nonlinear dynamics of bosonic junctions, described by the Gross-Pitaevskii equa- 
tion (GPE)£, is usually mapped to a simpler system characterized by two degrees of freedom 
(population imbalance and phase difference) while the nonlinear properties of the wave func- 
tion within the single well are neglected. In this approach the symmetric and antisymmetric 
stationary solutions of GPE are used as a basis to build a global wave function^. This 
description allows to show that for higher nonlinearities the symmetric solutions become 
unstable and degenerate to an asymmetric stationary (approximate) solution of the GPE 
corresponding to a new self-trapping regime^^. 

On the other hand, considering the double square well potential (instead of harmonic one), 
we discover that, in a wide range of nonlinearities, the system can either remain trapped 
mostly in one of the wells, or swing periodically from right to left and back. The switching 
from one state to the other is triggered by a slight local variation of the potential barrier 
between the wells. The coexistence of oscillatory and self-trapping regimes corresponds to 
the simultaneous presence of Josephson oscillations and of an asymmetric solution of the 
GPE. 

Our result differs from known behaviors of bosonic Josephson junctions, where the pres- 
ence of oscillatory or self-trapping regimes is uniquely determined by the parameters of 
the system. The resulting switching property should have a straightforward experimental 
realization in waveguide arrays, which constitute truly one- dimensional systems and are par- 
ticularly convenient for the observation of nonlinear effect j^ii&iiii&ii&iiLiS an d in engineered 
optical lattices of Bose-Einstein condensates^ 1 ^ 1 ^ 1 ^ 1 ^. 
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II. EXACT NONLINEAR SOLUTIONS IN DOUBLE SQUARE WELL 



Let us write GPE with double square well potential as follows: 

+ S " V(x)ij + H ^ = °' (1) 

where V(x) is the double square well (represented in Fig. [1]) with a total width 2L and the 
potential barrier height and width Vq and 21, respectively. 

The stationary solution of ([T]) are sought as ip(z, x) = Q(x) exp(— ij3z) with a real- valued 
function found in terms of Jacobi elliptic functions^ 

—L < x < -I : <E> = B cn[7 B (x + L) — K(k B ), k B ], 

I < x < L : $ = A cn[y A (x - L) + K(k A ), k A ], (2) 
—I < x <l : $ = C dn[7c(x - x ), k c ], 

with the parameters given in terms of the amplitudes by 

7a = v/^+A 7 B = v/^+A k\ = . A f t - , k 2 



c 2 

ll = w - (3 - k\ 



2 Vo-p-c 2 

2' c V -(3-C 2 /2-' 



where K denotes the complete elliptic integral of the first kind and by construction the above 
expressions verify the vanishing boundary values in x = ±L. 

The solutions are then given in terms of five parameters (A, B, C, Xq, f3), four of which 
are determined by the continuity conditions in x = ±Z. Thus the conserved total injected 
power (nonlinearity parameter) P t = J \ip\ 2 dx completely determines the solutions. Another 
useful conserved quantity is the total energy E given by 



E 



dip 



dx 



2 L/,14 



+ V(x)H 2 -— \dx. (3) 



In the weakly nonlinear limit (small Pt), the solutions are symmetric (odd or even). The 
even solution $+(x) corresponds to A = B in while odd solution corresponds to 

A = —B. For higher powers, namely above a threshold value, an asymmetric solution <& a (x) 
also exists for which A ^ ±B. These analytical solutions are represented in Fig. [TJ 

To plot the solutions we stick with the following parameter values: the width of the 
rectangular double well potential is 2L = 7.5, the barrier width is 21 = 0.25 and its height 
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FIG. 1: Plot of the double square well potential for the continuous model ([I]): 2L is the well width, 
Vo and 21 are barrier height and width. The curves are the plots of different types of solutions 
obtained for the total power Pt = 1.44. The inset shows the form of the asymmetric solution for 
different values of the total power. 

is Vo = 20. We derive the complete set of solutions ([2]) and display the dependence of their 
amplitudes on the total power P t = J \ijj\ 2 dx in the main plot of Fig. [2J Below the threshold 
value P t ~ 0.9 only the symmetric (odd and even) solutions exist and their amplitudes 
almost superpose. At the threshold value a new solution appears which is asymmetric with 
amplitudes A and B in the two wells, respectively, represented by the upper and lower 
branches in Fig. [2j 

III. TWO MODE APPROXIMATION 



The regime of Josephson oscillation is usually understood on the basis of coupled mode 
approach as follows. Using the symmetric and antisymmetric solutions, one builds a varia- 
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Total Power P 

FIG. 2: Dependence of the amplitudes (maximum values A and B of the expressions ([2])) of the 
symmetric and asymmetric solutions on the total power (the amplitudes of the odd and even 
symmetric solutions almost superpose). The inset displays the relative energy difference of the 
symmetric (<3?+) and assymetric (<3? a ) solutions in terms of the total power. 

tional anzatz by seeking the solution ijj(z,x) under the form 

ijj(z,x) = V>i(*)$i(z) + vfo{z)$2{x), (4) 
v/2$ 1 = $ + + $_, V2<&2 = <&+-<&_. 

The functions ^(z)! 2 and |?/; 2 (z)| 2 are interpreted as the probabilities to find the system 
localized either on the left or on the right part of the double square well. By construction, the 
overlap of $i with $2 is negligible, consequently, the projection of the GPE (OQ) successively 
on $1 and $2 provides the coupled mode equations^ 

oz 

+ D\^ 2 \ 2 ^ 2 = r^, (5) 
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with coupling constant r and nonlinearity parameter D defined by 



r = 



f[{d x $i){d x $ 2 ) + V$i$2] dx 
J ®l dx 



D 



J $\dx 
fQfdx' 



An explicit solution of (jSJ) in terms of Jacobi elliptic functions has been found ir*i and 
used in Bose-Einstein condensates iiA It is a good approximation for the system in a double 
harmonic potential welt^ and correctly describes the oscillatory regime in our case. Indeed, 
when the power is initially injected into one array, say |^i(0)| = 1, ^(O)] = 0, we obtain 
for D < 4r 



Since oscillates around the value 0, this expression describes an oscillation of light 
intensity between the left and the right wells. The period of this oscillation is 



and has been checked on various numerical shots at different total input power. In summary, 
while the self-trapping regime is directly interpreted in terms of the asymmetric solution, 
the interpretation of the Josephson oscillation regime needs to call to the coupled mode 
approach, which in turn fails to explain the observed coexistence of both regimes. 

Such a coexistence, however, is understood in terms of the energy ([3]) which can be eval- 
uated, at given total power P t , both for the symmetric solution $ + and for the asymmetric 
solution $ a . As shown in the inset of Fig. [2] these two energy values E + and E a turn out to 
be very close up to the total power value Pt ~ 2. Consequently, switching from a regime to 
the other is allowed at fixed power. In particular, in the numerical experiments of Fig. [3l 
total power and energy are the same before and after the local variation of the potential 
barrier value. 

It is worth to remark that a similar analysis in the case of harmonic double well 
potentials'^ shows that the energy of the asymmetric solution (when this solution ex- 
ists) is significantly smaller than the energy of the symmetric solution. In such a situation, 
it is thus impossible to switch from a self-trapped state to an oscillatory regime when keeping 
both the energy and the total power constant. 




(6) 



T = 2K(D/4r)/r 



(7) 
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FIG. 3: Numerical simulation of the GPE equation ([I]). By a slight local variation at z = 150 
of the potential barrier height, represented in the inset, the regime switches from self-trapping to 
Josephson oscillations. The injected total power is Pt = ^ IV'jl 2 = 1-44. 

IV. APPLICATIONS FOR BEC AND COUPLED WAVEGUIDE ARRAYS 

Now our aim in this section is to suggest the realistic experiments on BEC and waveguide 
arrays, which are engineered in such a way to mimic two weakly coupled chains of JJ's (see 
Fig. H]). Using such an experimental set-up, we demonstrate the feasibility of the efficient 
control of a switch between oscillating and self- trapping states of the systems. We show 
that our problem reduces to the Gross-Pitaevskii equation (GPE)^ in a double square well, 
which displays very different properties from the previously considered double harmonic well 
potential^ 1 ^^^ 1 ^. Our results are broadly applicable and open the way to the experimental 
study of these phenomena in the dynamics of those systems. 
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FIG. 4: Schematics of the suggested experimental setup. In the context of BEC the optical lattice is 
supplemented by large energy barriers from both sides and a small one in the middle (solid curve) . 
The condensate is initially loaded mainly into the right part of the optical lattice (the dashed line 
represents particle density). The inset shows the reduction of the problem to the particle motion 
in a double square well potential (details are given in the text). In the context of waveguide arrays 
the solid curve displays (with the opposite sign) variation of the refractive index across the array. 

We start the consideration of the case of BEC in an optical lattice, for which a one- 
dimensional Hamiltonian has a following form: 

dib h 2 d 2 ib . . 2h 2 a s . , , 9 , 

where m is atomic mass, a s < is the scattering length corresponding to the attractive atom- 
atom interactions and a± = \Jh/muj±_ is the transversal oscillation length, which implicitly 
takes into account the real three dimensionality of the system^, u± being the transversal 
frequency of the trap. The optical lattice potential is 

V(x) = v cos 2 (k L x) for \ki,x\ > n/2 
V(x) = (v + V ) cos 2 (k L x) for \k L x\ < tt/2, (9) 

where ki is the wavenumber of the laser beams that create the optical lattice and Vq is the 
height of the additional spatial energy barrier placed in the middle of the optical lattice. 
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v = — > = 9 = — 4 , (11) 



Besides that, Dirichlet boundary conditions with if}(±L) = are chosen in order to describe 
the large confining barriers at both ends of the BEC. These boundary conditions could be 
realized experimentally by an additional optical lattice with larger amplitude and larger 
lattice constant, as shown in FigJH 

Introducing a dimensionless length scale x = 2k B x and time t = Egt/h, where Eb = 
8Er = Ah 2 k 2 L /m and Er is the recoil energy^, we can rewrite (Tl8|) as follows 

ld 2 ^ 

i ^ = —2-d^ + n ^ + mH! ' (10) 
where the normalized wave-function, j \ty(x)\ 2 dx = 1, is introduced^. The dimensionless 
potential V still has the form ([9]) with the following dimensionless depths of the optical 
lattice 

v ~ _ V _ Na s 

E B E B k h a\ 

g being the dimensionless nonlinearity parameter. 

We have performed numerical simulations of Eq. (flUj) with 12 wells (6 wells on each side 
of the barrier as presented in Fig. H]) and the parameters v = 0.25 (in physical units this 
means that the depth of the optical lattice is v = 2Efj), Vq = 0.15 and we fix the nonlinearity 
to the value g = —0.025, i.e., we choose attractive interactions. The dynamics is similar for 
repulsive interatomic forces (see the discussion below). The phenomenon we study in this 
Letter does not depend significantly on the actual size of the system, if at least 3 lattice 
sites are present at each side of the barrier. 

As seen from the left panel of Fig. [5j if one prepares the condensate in a self-trapped state 
it remains there until we apply the pulselike time variation of the barrier displayed in the 
inset. After that action, the condensate goes into the oscillating tunneling regime. On the 
other hand, preparing the condensate in the oscillating tunneling regime (right graph in Fig. 
|5]) one can easily arrive at a self-trapping state by varying again the energy barrier in the 
middle as displayed in the inset. Let us mention that, as far as the energy of the barrier is 
changed adiabatically, the total energy of the condensate does not vary, i.e. the self-trapped 
and tunneling oscillatory regimes have the same energy. This is quite different from what 
happens in a double harmonic well potentia^ 1 ^*^ 1 ^' 1 ^. The point is that, in the double 
harmonic well, the asymmetric stationary solution is characterized by a smaller energy than 
the symmetric solution and this difference increases sharply with increasing nonlinearity. 
Hence, a drastic energy injection is required in order to realize the transition between the 
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FIG. 5: Numerical simulations of Eq. (|1U|) : the left graph represents the transition from a self- 
trapped state to the macroscopic tunneling regime, while the right graph describes the inverse 
process. The insets in both graphs show the variation of the energy barrier necessary to realize the 
switching between the different regimes. 

two regimes; whilst in our case the transition is simply achieved only by varying pulsewise 
the energy barrier. Below we argue that this happens because our case effectively reduces 
to the case of a double square well potential (see the inset of Fig. H] and the reduction 
procedure below) for which asymmetric and symmetric stationary solutions carry almost 
the same energies in a wide range of the nonlinearity parameter. 

Now we proceed to reducing Eq. (jTUl) to a Discrete NonLinear Schrodinger equation 
(DNLS). We discretize it via a tight-binding approximation^^^, representing the wave 
function ty(x) as 



where <&j(x) is a normalized isolated wave function in an optical lattice in the fully linear 
case g = and could be expressed in terms of Wannier functions (see, e.g.,—). For clarity, we 
use here its approximation for a harmonic trap centered at the points rj = jir(\j\ + l/2)/\j\ 
varies from 1 to n, the number of wells). In the context of the evolution equation (llOjl 
<&j(x) has the form 



for \j\ 1, and one should substitute v by v + Vq in the above expression in order to get an 




(12) 



j 




(13) 
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approximate formula for the wave function for \j\ = 1. 

Assuming further that the overlap of the wave functions in neighboring sites is small, we 
get from (TTUI) the following DNLS equation for the sites \j\ ^ 1 



dt 



(14) 



while for 



1 we have 



ill 



d<f>.. 



dt 



-Q4>±2 - Qi<j> T i + Ui\<p±i\ 2 <f)±i, 



(15) 



where we assume pinned boundary conditions. The constants Q, Qi, U and U\ are easily 
computed from the following expressions (|j| 7^ 0): 



Q = - 
Qi = - 



dx dx 
cte dx 



cos (x/2)$ i $ J+ i 



+ (0 + V r )cos^(x/2)$ 1 <1) 



dx, 



U = g J dx ~ Ui = g J $± 2 dx 



(16) 



In order to characterize the solutions of Eqs. ( lT4l) and ( 1T5|) . we follow the same procedure 
used in Ref.— , which goes through a continuum approximation. Assuming that 0i = 0_i 
we finally arrive at 

ihd<f>{j) d 2 (p(j) 



+ W(j)cj>(j) + R\(f>(j)\*(j>(j), 



(17) 



where now j is a continuous variable, is a double square well potential with a barrier 

height w = 2{Q — Qi)/Q and width / = 1, <p(j) obeys pinned boundary conditions <p{j = 
±L) = (2L is a width of a double square well potential) and the nonlinearity parameter 
is given by R — U/Q < 0. Expressing -0 = \J\R\(j) and 2 = Qt/^ and mentioning that total 
power P t is connected with nonlinearity parameter R as P t = \R\, we see the equation (TTTf) 
is the same as ([1]) and thus all the above consideration of peculiarities of double square well 
potential directly applies to the considered BEC lattices. 

In case of the waveguide systems the situation is even simpler. Particularly, as well known 
an array of adjacent waveguides coupled by power exchange is modeled by the discrete 
nonlinear Schrodinger equation (DNLS)^ 1 ^ which reads 



+ —{rij - n)ipj + Q(ip j+1 + - 2ipj) + \ipj\ 2 ipj 



dz 



0. 



(18) 
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where waveguides discrete positions are labelled by the index j (-N < j < N), and the 
complex field results from the projection of the electric field envelope on the eigenmode of 
the individual waveguide. It is normalized to a unit onsite nonlinearity. The linear refractive 
index rij is set to n for all j ^ 0, and to n < n for j = 0. The coupling constant between 
two adjacent waveguides is Q and u and c are the light frequency and velocity. Vanishing 
boundary conditions ipN+i = 4>-n-\ = model a strongly evanescent field outside the 
waveguides. Considering now 1 / ^[Q as being a virtual grid spacing we may represent ipj (z) 
by the function ip(x, z) in the continuous variable x = j / \[Q. As a result the DNLS model 
( TT8l) maps to the (CQ) with a double square well potential considered initially. 

V. CONCLUSIONS 

A new coherent state in square double well potential has been discovered. This coherent 
state has the property of being bistable: one can easily switch from oscillatory to self- 
trapping regimes and back. This nontrivial behavior may have interesting applications in 
various weakly linked extended systems, such as Bose-Einstein condensates, waveguide or 
Josephson junctions arrays, which deserve further studies. 

In the region of nonlinearities where the asymmetric solution coexists with the symmetric 
and asymmetric stationary solutions, we have induced the switch from one regime to the 
other by varying the height of the barrier. In view of a real experiment one could induce 
such flips by varying the refractive index of the central waveguide (in the context of weakly 
linked waveguide arrays) or by pulswize change of optical barrier potential (in case of BEC). 
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